Path-integral evolution of chaos embedded in noise: Duffing neocortical analog 
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Abstract — A two dimensional time-dependent Duffing oscillator model of macroscopic neocortex 
exhibits chaos for some ranges of parameters. We embed this model in moderate noise, typical of the 
context presented in real neocortex, using PATHINT, a non-Monte-Carlo path-integral algorithm that is 
particularly adept in handling nonlinear Fokker-Planck systems. This approach shows promise to 
investigate whether chaos in neocortex, as predicted by such models, can survive in noisy contexts. 
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1. INTRODUCTION 

A global theory of neocortical dynamic function at macroscopic scales was developed to explain 
sahent features of electroencephalographic(EEG) data recorded from human scalp [1]. In some limiting 
cases, the theory predicts EEG standing waves composed of linear combinations of spatial modes. These 
modes (which may be called order parameters in the parlance of modem dynamical methods) are 
governed by Unear or quasilinear ordinary differential equations [1-3]. These data and theoretical works 
suggest several connections to reports of chaotic attractors estimated for EEG data [1]. In order to 
illuminate possible relations between sub-macroscopic chaos (for example in neocortical columns) and 
large scale data observed at the scalp, we suggested a simple metaphoric system, the linear stretched 
string with nonlinear attached springs [4-6]. We have studied a forced partial differential equation 
describing the string-spring system, where the forcing term might represent a steady-state evoked 
potential in the neocortical analogy. In the limiting case when string tension is zero, the forced Duffing 
equation is obtained. The Duffing equation exhibits a rich dynamic behavior of periodic and chaotic 
motion, for different parameters (associated with different wave numbers in the string-spring system) [7]. 

Another series of papers has outHned a statistical mechanics of neocortical interactions (SMNI), 
demonstrating that the statistics of synaptic and neuronal interactions develops a stochastic background of 
moderate noise that survives up to the macroscopic scales that are typically measured by scalp 
EEG [4,8-12]. 

This study is concerned with the relation between chaos that may occur in particular spatial modes 
(the order parameters) and EEG signals consisting of combinations of modes and noise. In this context, 
we may consider the "noise" to be either instrumental or biologic, e.g., as detailed in the SMNI papers. 
That is, can we reasonably expect correlation dimension estimates of EEG data to provide information 
that can be used to construct physiologically-based theories of neocortical dynamics. 

Section 2 describes the Duffing metaphor of neocortex that exhibits chaos. Section 3 describes the 
embedding of the deterministic model into a stochastic medium. Section 4 describes the PATHINT code 
used for the stochastic model. Section 5 gives the results of the PATHINT calculations. Section 6 is a 
brief conclusion. 



Path-integral of Duffing model 



-3- 



Ingber, Srinivasan and Nunez 



2. DUFFING ANALOG OF NEOCORTEX 

At this stage of our understanding, choosing a specific fit to a macroscopic EEG model is 
necessarily limited to the specific data and would necessarily be controversial. However, a relatively 
simple physical system in which the general issues of spatial scale interaction and connection to 
experimental EEG can be discussed with minimal complication is a linear stretched string with attached 
nonlinear springs [4]. 

The ends of the string may be considered fixed or the string may form a closed loop to make the 
system more analogous to the closed neocortical/corticocortical fiber surface. A closed loop of 
transmission fine with nonlinear dielectric conductance is the equivalent electric system. In the proposed 
metaphorical Hnk with the brain, the springs represent local circuit dynamics at columnar scales whereas 
the string is beUeved to be somewhat more analogous to neocortical interactions along the corticocortical 
fibers [5,10]. 

String displacement, considered to be somewhat analogous to the modulation of cortical surface 
potential is governed by 

In the case of fixed-end supports, the forcing function takes the form 

Fix, t) = B(x) cos t . (2) 

Here is the usual string tension parameter (or wave velocity parameter), Wq is proportional to the linear 
spring constant, and a is the damping parameter. In the neocortical metaphor, F(x, t) might represent 
ongoing random subcortical input (e.g., from the thalamus), or sensory information designed to study 
Hnear/nonlinear resonance. Since the string/springs is simply a metaphor for the brain, the nonhnearity is 
arbitrary, so we choose a cubic nonhnearity and sinusoidal driving so that string displacement is governed 
by the partial differential equation 

- ^ V ^ " 37 ^ ^ ^ = ^^'^ '''' ' ■ 

The driving has been included in order to connect the system to the well studied forced Buffing's equation 
as well as to obtain chaotic dynamics. In the hmiting case of a completely relaxed string (i.e., the velocity 
of propagation is zero), the springs behave independently and drive the string, governed by the forced 
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Buffing's equation atx = xq. 



+ a ^ + ci)Q<^ + <^= B(x) cos t . 



(4) 



The alternative way to arrive at Buffing's equation is to do a spatial mode expansion of the string equation 
and then only keep the lowest wavenumber [6]. For instance, the Lorenz system is derived in this way 
from the atmospheric turbulence equations [13]. 

This equation has been studied numerically for the case a)Q = in the parameter space (a, B) using 
spectral analysis and phase -plane projections [7]. A rich structure of periodic, quasi -periodic and chaotic 
motion were reported. We duplicated Ueda's results using a standard fifth and sixth order Runge-Kutta 
algorithm (IMSL routine BIVPRK). The dynamics were classified as chaotic on the basis of 
characteristic broadening of spectral lines and aperiodic trajectories in phase space as well as correlation 
dimension analysis. In all simulations discussed here, the damping parameter a = 0. 05 and B = 5.5. 

Figs, la and lb illustrate the evolution of the phase space for the non-chaos and chaos cases of 
»o = 1-0 and (Oq = Q.\ cases, resp., after transients have been passed >250). Figs. 2a and 2b illustrate 
the evolution of the phase space for the non-chaos and chaos cases of »o = 1- and o}q = 0.\ cases, resp., 
during a transient period between t = l2 and f = 15.5. The non-chaos cases are similar in appearance 
except for the settling of transients. 



3. EMBEDDING DETERMINISTIC DUFFING SYSTEM IN NOISE 



The deterministic Buffing system described above can be written in the form 



X = fix, t) 



f = -ax - WqX + B cos t 



(5) 



and recast as 



x = y 



y = fix, t) , 

f = -ay - coqX + B cos t . (6) 
Note that this is equivalent to a 3-dimensional autonomous set of equations, e.g., replacing cos / by cos z. 
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defining z = P, and taking p to be some constant. 

We now add independent Gaussian-Markovian ("white") noise to both x and y, rj\, where the 
variables are represented by i = {x, y} and the noise terms are represented by j = {l, 2}, 

x = y + g]crji , 

y = fix, t) + 15772 , 

< r]\t) >,= , 

<rj\t),rj^\t')>r,= 5^^'S{t-t') . (7) 

In this study, we take moderate noise and simply set gj = l. OS- . 

The above defines a set of Langevin rate-equations, which can be recast into equivalent Fokker- 
Planck and path-integral equations [14]. A compact derivation has been given in several papers [10,15]. 
The equivalent short-time conditional probabihty distribution P, in terms of its Lagrangian L, 
corresponding to these Langevin rate-equations is 

1 

P[x, y; t + At\x, y, t] = ^ exp(-LAf) , 



4. PATfflNTCODE 

PATHINT is a non-Monte-Carlo histogram C-code developed to evolve an n-dimensional system 
(subject to machine constraints), based on a generahzation of an algorithm demonstrated by Wehner and 
Wolfer to be extremely robust for nonlinear Lagrangians with moderate noise [16-18]. This algorithm has 
been used for studies in neuroscience [19,20], financial markets [21,22], combat analysis [23,24], and in a 
study of stochastic chaos [25]. 

PATHINT computes the path-integral of a system in terms of its Lagrangian L, 

( t \ 



S{q{tQ) = qo)Siqit) = q^) , 



to J 
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u+l 

Dq = lim n g''^ n (2^Atr"^dq'p , 



L{q',q',t) = ^-{q' - g')gu<q' - g') , 



gii (§ ) ' 

g = det(g,r) • (9) 

If the diffusions are not constant, then there are additional terms [14]. 

The histogram procedure recognizes that the distribution can be numerically approximated to a high 
degree of accuracy as sums of rectangles at points q^ of height P, and width Aq'. For convenience, just 
consider a one-dimensional system. The above path-integral representation can be rewritten, for each of 
its intermediate integrals, as 



Fix; t + At) = j dx'[g^''^{2jtAtr^'^ exp(-L At)\P{x; t) 

= J dx'G(x, x'; At)P{x'; t) , 



P{x; t) = j:^(x-x')Pi{t), 

i=l 



7r(x-x') = • 



1 , ix' - - Ax'-^) <x<ix' + - Ax') , 

2 2 (10) 

0, otherwise. 



This yields 

Pi{t + At) = Ty(At)Pjit), 

TJAt) = r f ^"+A^'/2 dx { xUAxJ/2 dxG{x, X\ At) . (1 1) 

■' Ax'-^ + Ax' J x'-Ax'-yi J xi-Axi-y2 

Tjj is a banded matrix representing the Gaussian nature of the short-time probability centered about the 
(possibly time-dependent) drift. 

Here, this histogram procedure s used in two dimensions, i.e., using a matrix Tfj/^i [23]. Explicit 
dependence of L on time t also can be included without complications. Care must be used in developing 
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the mesh in A^', which is strongly dependent on the diagonal elements of the diffusion matrix, e.g., 

^q' = {^tg''f'^ . (12) 

Presently, this constrains the dependence of the covariance of each variable to be a (nonlinear) function of 
that variable, in order to present a straightforward rectangular underlying mesh. 

Since integration is inherently a smoothing process [21], fitting data with the short-time probability 
distribution, effectively using an integral over this epoch, permits the use of coarser meshes than the 
corresponding stochastic differential equation. For example, the coarser resolution is appropriate, 
typically required, for numerical solution of the time-dependent path integral. By considering the 
contributions to the first and second moments conditions on the time, variable meshes can be derived [16]. 
The time slice essentially is determined by A?<L~', where L is the uniform Lagrangian, respecting 
ranges giving the most important contributions to the probability distribution P. Thus, A? is roughly 
measured by the diffusion divided by the square of the drift. 

Monte Carlo algorithms for path integrals are well known to have extreme difficulty in evolving 
nonlinear systems with multiple optima [26], but this algorithm does very well on such systems. The 
PATHINT code was tested against the test problems given in previous one-dimensional systems [16,17]. 
Two-dimensional runs were tested by using cross products of one-dimensional examples whose analytic 
solutions are known. 

5. PATfflNT CALCULATION OF DUFFING SYSTEM 

While we appreciate that a true calculation of chaos of a stochastic system is quite difficult, here 
our approach is to simply perform a rigorous numerical calculation of the evolution of the Fokker-Plank 
equation in its equivalent path-integral representation. Then, we can compare the chaos and non-chaos 
topologies of the deterministic system, which have been seen above to (barely) be different, to the 
topologies of the stochastic system. The topology of chaotic systems is an important invariant that often 
can be used to identify its presence [27]. 

A mesh of A? = 0. 1 within ranges of {x, y} of ±25 was reasonable to calculate the evolution of this 
system for all runs. To be sure of accuracy in the calculations, off-diagonal spreads of meshes were taken 
as ±5 increments. This lead to an initial four-dimensional matrix of 159 X 159 xllxll = 3 059 001 
points, which was cut down to a kernel of 2 954 961 points because the off-diagonal points did not cross 



Path-integral of Duffing model 



-8- 



Ingber, Srinivasan and Nunez 



the boundaries. Reflecting Neumann boundary conditions were imposed by the method of images, 
consisting of a point image plus a continuous set of images leading to an error function [28], but 
calculations support the premise that they are unimportant here as the evolving distributions stayed well 
within the boundaries. Since we have added noise to a set of ordinary differential equations, bringing 
along additional boundary conditions, we thought it prudent to establish solutions that would be rather 
independent of any Dirichlet or Neumann boundary conditions. 

A Convex 120 supercomputer was used, but there were problems with its C compiler, so gcc 
version 2.60 was built and used. Runs across several machines, e.g.. Suns, Dec workstations, and Grays, 
checked reproducibility of this compiler on this problem. It required about 12 CPU min to build the 
kernel and perform the calculation of the next distribution at each A?-folding. Data was accumulated after 
every 5 foldings. The chaos case was run up until t = 33, requiring 67.5 CPU hours, and the non-chaos 
case was run up until t = 15. 5, requiring 32.0 CPU hours. 

The above calculations took over 100 CPU-hours, and this was a self-imposed hmit for this first 
examination of this problem. As demonstrated by the deterministic calculations above, further 
calculations would have to be performed to enter into the time domain of ? > 50 after transients have 
diminished and where chaos truly is strongly exhibited. Furthermore, for t > 20, we noticed that there is 
sufficient diffusion of the system towards the boundaries at {x, >'} = ± 25, so that these ranges would have 
to be increased, thereby consuming even more CPU time per At run. For the purposes of this paper, we 
just look at epochs where the deterministic calculations show differences between the chaos and non- 
chaos cases, and compare these to our observations of the stochastic system. 

Figs. 3a, 3b and 3c illustrate the stochastic evolution of the phase space for the non-chaos case of 
©0 = 1-0 during a transient period of t=l2- 15; compare to Fig. 2a. Figs. 4a, 4b and 4c illustrate the 
stochastic evolution of the phase space for the chaos case of »o = 0- 1 during a transient period of 
? = 12- 15; compare to Fig. 2b. The differences in the stochastic cases are practically non-existent, and 
any differences that are present are likely due to differences in drifts due to changes in the parameters. 

6. CONCLUSION 

We have presented calculations of a two-dimensional time-dependent Duffing system, that has some 
support for modeling interactions in macroscopic neocortex, that possesses chaotic and non-chaotic 
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regions of activity. We embedded this system in moderate noise and investigated the epochs within the 
transient periods of these systems to see if there are any apparent differences between the deterministic 
and stochastic systems. 

Future studies taking hundreds of CPU-hours will investigate time epochs beyond the transients to 
see if moderate noise suppresses chaos in such models of neocortex. 
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FIGURE CAPTIONS 

FIG. 1 . Phase space plots of <l> versus ^. (a) The evolution of the phase-space for the non-chaos 
case of coq = 1.0, after transients have been passed, (b) The evolution of the phase-space for the chaos 
case of (Oq = 0A, after transients have been passed. 

FIG. 2. Phase space plots of <I> versus ^. (a) The evolution of the phase-space for the non-chaos 
case of ©0 = 1 • 0' during the transient period of 12 - 15. 5. (b) The evolution of the phase- space for the 
chaos case of (Oq = 0A, during the transient period of?=12-15.5. 

FIG. 3. Conditional probability distributions as a function of <I> versus <i>. (a) The stochastic 
evolution of the phase space for the non-chaos case of »o = 1 • at ? = 12. (b) The stochastic evolution of 
the phase space for the non-chaos case of t»o = 1- ^ = 13. 5. (c) The stochastic evolution of the phase 
space for the non-chaos case of «o = 1- at ^ = 15. 

FIG. 4. Conditional probability distributions as a function of <J> versus (a) The stochastic 
evolution of the phase space for the chaos case of a)Q = 0.l at t = l2. (b) The stochastic evolution of the 
phase space for the chaos case of (^o = 0. 1 at ? = 13. 5. (c) The stochastic evolution of the phase space for 
the chaos case of = 0. 1 at f = 15. 
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Figure la 
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Figure lb 
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Figure 2a 
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Figure 2b 
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Figure 3a 
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Figure 3b 
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Figure 3c 
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Figure 4a 
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Figure 4b 
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Figure 4c 




